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An accurate and fast sampling method was developed on the modeling of a voxel phantom called Rad- 
HUMAN for radiation protection of MC-based radiation transport and simulation. The segmented organ voxels, 
which were assigned with three dimensional (3D) coordinates, were simplified through a two-step hybrid sam- 
pling algorithm. Firstly, certain voxels were sampled into a coordinate matrix by the nearest neighbor sampling. 
Secondly, the coordinate matrix was renewed using a weighted sampling. To compare visualization with the 
sampling, the resultant matrix was used to extract the contour of organs/tissues for constructing polygon-surface 
phantom. The feasibility and effectiveness of the sampling method was verified through the modeling of large 
organs (e.g. skeleton system) and application of transformation to MC computational geometries. 
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I. INTRODUCTION 


Computational phantoms have been used extensively in 
radiation protection for Monte Carlo (MC)-based radiation 
transport and simulation [1]. Accuracy of the phantom plays 
a critical role in dose calculation [2, 3], and computational 
phantoms based on color photographic images and contain- 
ing detailed information of human body have been developed 
in many countries [4, 6]. However, it is a big problem to deal 
with the massive voxels accurately [5]. Lots of refinement 
methods were developed not directly focusing on the origi- 
nal voxels, and sometimes voxels were sampled with inter- 
polation. This may generate some ambiguity that may cause 
serious problem in whole-body phantom construction. Due 
to limitation of computer performance, the large organs (e.g. 
skeleton system) from the image slices cannot be constructed 
as a whole accurately. 

Based on high-resolution sectioned images of a Chi- 
nese Visible Human (CVH) dataset, a voxel phantom called 
Rad-HUMAN (Accurate whole-body computational phan- 
tom of Chinese adult female) was established by the FDS 
Team (www. fds.org.cn). 

This CVH data set, which contains 3641 slices with 
3286 x 1586 pixel resolution, was obtained from a Chi- 
nese female cadaver. Each voxel size of the data set is 
0.15mm x 0.15mm x 0.25mm for the head and neck and 
0.15mm x 0.15mm x 0.5mm for the rest of the body. The 
total number of the voxels is about 16.8 billion. Segmen- 
tation [6] (Fig. 1) was processed manually and calibrated 


* Supported by National Special Program for ITER (No.2011GB113006), 
Strategic Priority Research Program of Chinese Academy of Sciences 
(No.XDA03040000), National Natural Science Foundation of China (No. 
91026004), the Knowledge Innovation Projects of Chinese Academy of 
Sciences (No.095CF2R211) and the Key Foundation for Young Talents in 
College of Anhui Province (No.2013SQRL063ZD) 

t Corresponding author, yican.wu @fds.org.cn 


by MCAM (Multi-Physics Coupling Analysis Modeling Pro- 
gram) developed by the FDS Team [7-9]. In addition to the 
fact that the voxel phantom was difficult for Monte Carlo 
dose calculation, demonstrating its geometry was also labo- 
rious [6]. 


Fig. 1. (Color online) Organ segmentation. 


In this work, a flexible sampling method was investigated 
on the modeling of Rad-HUMAN voxel phantom. Based on 
the method, the visualization was compared and computa- 
tional geometries for Monte Carlo dose calculation was dis- 
cussed. 


Il. VOXEL SAMPLING METHOD 
A. Nearest neighbor sampling 


After segmentation, the organs had different RGB colors to 
be distinguished. In this paper, the data set of voxels assigned 
with 3D coordinates can be written as 


(0,0, h — 1) (w — 1,0,h — 1) 


(dij k)wxixh = : : 
(0, -—1,h—1) (w—1,l—1,h—-1) 
i,j € [0,w — 1/1 — 1); 
(1) 
where w, l and h are the width, length and height of the CVH 


dataset, respectively; and 2, j and k are the position subscripts 
of the voxels. 
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It is no use obtaining all the voxels of a certain organ to 
arduously construct its model. In the nearest neighbor sam- 
pling, the segmented voxels are appropriately enlarged and 
sampled from their nearest neighbors. 

Change the voxel array from w x 1 x n to w’ x I! x n’, the 
sampling enlargement factor can be expressed as 


w l n 

(Sz, Sy, Sz) — (or? J nh: (2) 

To be in consistence with the later sampling, the voxel co- 
ordinate of an organ can be described as a sparse matrix a 


a= (ai jk)wxixh, 
pw Í digem SF (diga) K 
Mi = 1 0, BF (dijk), 


where AF (di j,k \(OrÅF (di j k)) means that d; j,ẹ is (or not) a 
coordinate of the original organ. 

In the nearest neighbor sampling, the sampled organ 
f(a, y, Z) is acquired from an original organ F(X,Y, Z) de- 
scribed as 


f(x,y, z) 


where symbol “[ ]” stands for rounding off the number for 
the result. So the resultant coordinate matrix of the sampled 
organ can be described as 


= F([Se x 2], [Sy x y] [52x A) &) 


Y = (ijk Jwr xt xh! 
ijk = di j,k, J F(duxs, },[9x Sy], [kx Sz ) (5) 
ij, 0; Aldiss tesco eat): 


To decrease the demand for computer memory, in present 
study, the nonzero coordinate (Yi j,k) was stored on disk as 
file format corresponded with exclusive address for com- 
pressed storage. 


B. Weighted sampling 


As accuracy loss is serious by the nearest sampling, a 
weighted sampling in the present study was proposed. With 
the same enlargement factor defined in Eq. (2), a sampling 
unit can be enclosed by a lattice with a size of Sy x Sy x Sy. 
Obviously the voxels in the lattice cannot always belong to 
one organ, hence a weighted factor to judge whether the vox- 
els of the sampling unit can be regard as a voxel after sam- 
pling. 

After sampling, the relative position of the voxel was 
changed corresponding with the sampled voxel array (w’ x 
U x n’). That is, the coordinates must be zoomed by a size 
of (1/S,) x (1/5,) x (1/S,) for the consistence. In addi- 
tion, described as the coordinates with exclusive addresses, 
the original voxels of the lattice can be projected to a new ad- 
dress of the sampled storage. This brought the effectiveness 
for numbering the amounts of the voxels which belong to one 
organ being projected. So according to a defined weight of 
the voxels in the lattice, the voxels of a certain organ can be 
sampled fast. The detailed steps for the weighted sampling 
are as follows. 
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Step 1: Obtain coordinate matrix a of an organ referred as 


Eq. (3). 


Step 2: Map the matrix a to a sampling matrix 3. Supposing 


(La; Yas Za) is an element of a, the element of 8 be- 


comes 
(xB, YB, ZB) E (LasYasZa) W 
1/S, 0 0 
o=| 0 Is, 0. h (6) 
0 oO 1S, 


where the mapped address in 6 can be described as 0 
with the element obtained from 


0 = (Las Ya; Za) TT 
w wxl ) (7) 


Sy? SaxX8y? Sy X8y XS 


T = 


w and 7 in Eq. (7) are factors to multiply every element 
of a. 


Step 3: Record the weights of certain coordinates which are 


mapped into the same locations. 
Let (x, y, z) be a coordinate after sampling, which are 
mapped from o (a region of sampling lattice). 
a €{([x x Sz], [y x Sy], [z x Sz]), 
(Mæ +1) x Se], [y+ 1) x Sy], (8) 
[(2 + 1) x Sz})}- 


In this step, the coordinates in o is stored at the same 
location with (a, y, z) recording their amounts. 


Step 4: Filter the coordinates from 8, with the weights of no 


less than a weight factor. 


Supposing that the weights were defined as the propor- 
tional amounts of the sampling unit which is 0.5, the 
filtered coordinate matrix y can be selected from 8 with 
the condition of 


(£, y, zZ) > 0.5 x Sz x Sy x Sz. (9) 


The resultant coordinate matrix y can also be stored on 
disk as file format in sequential order. 


C. Hybrid voxel sampling 


The nearest sampling is a fast way to reduce the data size. 
However, it is not recommended to use this sampling method 
on some smaller data because of serious accuracy losses. In 
this paper, we combine it with the weighted sampling as a 
two-step hybrid sampling algorithm. In the hybrid sampling, 
the sampling unit is divided into smaller lattices for the near- 
est sampling, and then the smaller lattice can be assighed with 
weights for the weighted sampling. 

To that end, the enlargement factor of (Sz, Sy, Sz) was 
divided into (Nz, Ny, Nz) and (P,, Py, P2) for the nearest 
sampling and the proportional sampling, respectively. 
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(Po, Py, Pa) = (F Sy =) (10) 


Nr’ Ny’ Nz 
Using an enlargement factor of (Nz, Ny, Nz) for the near- 
est sampling, the original coordinate matrix can be obtained 
according to Eq. (5) where 


di j,k, SE (dex Na], [jx Ny] [kx No]) 
ijk = H a un 11 
Oe { 0, AF (dux Nel ix Nyl,[&x N1}: eh 


In order to combine the nearest sampling with the weighted 
sampling for the hybrid sampling, the factor in Eqs. (6) 
and (7) for the mapping and locating becomes 


1/P, 0 0 
w=[{ 0 1/P, 0 (12) 
0 0 1/R 


and T = (1/Pp, w/(P} X Sz), w x l/(P; X Sz X Sy)), re- 
spectively. 

The detailed steps to implement the hybrid sampling algo- 
rithm are shown in Fig. 2 . 


(Begin sampling 


(Px,Py,Pz) <3 (Sx/Nx,Sy/Ny,Sz/Nz) 


| 


Obtain a by nearest sampling using (Nx, Ny, Nz) 


| 


B (3a: ©, with (x,y,z) recorded | 


| 


ly p, with 5(x,y,z)>0.5xPxxPyxPz | 


~ = 


oe 
(End sampling ) 


Fig. 2. Hybrid voxel sampling algorithm. 


Il. RESULTS AND DISCUSSION 
A. Application in construction of polygon-surface phantom 


The organs or tissues of interest (e.g., lungs, liver, skin 
etc.) from the original tomographic photography were identi- 
fied by assigning every pixel with an identification numbers. 
All these numbers can be stored sequentially as a pixel ma- 
trix that can be used to extract equivalent matrix, from which 
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TABLE 1. The comparison of the heart reconstruction 


Comparison Nearest Weighted 
Enlargement factor proportional 2.0 2.0 2.0 2.0 2.0 2.0 
Weight 0.5 

Voxel size (mm?) 0.3 x 0.3 x 1.0 0.3 x 0.3 x 1.0 
Times (s) 19 268 

Number of triangular faces 589 005 971133 

Holes 1516 1274 

Model size (MB) 29.2 28.5 


the polygon-surface model can be constructed by Marching 
Cubes [10]. 

The equivalent matrix can be acquired by the last step of 
the sampling method. For example in weighted sampling, on 
purpose of generating the matrix y, the coordinates satisfying 
the condition of d(x, y,z) > 0.5 x Ss x Sy x Sz are set as 
an identification number, while the coordinates not satisfying 
that condition are set as another identification number. The 
changed matrix y can then be used as the equivalent matrix 
using VTK toolkit to generate the polygon-surface model. 

Preparatorily, these sampling methods were processed with 
the just one CPU of the 3.10 GHz Intel Core™ 2Quad Pro- 
cessor 15-2400 64 bit operating system. To compare with two 
special cases of the hybrid sampling, the construction of heart 
organ containing 265 slices was taken for example. Defin- 
ing the proportional amounts as the weights, Table | presents 
relevant sampling parameters between the nearest and the 
weighted sampling methods which are compared by the vi- 
sualization in Figs. 3(a) and 3(b). 

When the voxel increased with a factor of 2.0 in x, y and z 
directions, the time expenditure by the nearest sampling was 
about 7.09% of that by the proportional sampling. The num- 
ber of the holes and model size produced by the weighted 
sampling was 84% and 95% of the one produced by the near- 
est sampling, respectively. By the weighted sampling, the 
incorrect model incurred by segmentation can be repaired 
(Fig. 3(d)). 

Parameters between the two methods (Table 1 ) are mainly 
dependent on computer performance and data set resolution 
etc. For precision adjustment, the proportional weight factor 
can be used. 


B. Improvement on the sampling and analysis 


The nearest and the weighted sampling method in present 
paper are implied as two special cases of the hybrid sam- 
pling. To discuss their combined criteria, extra experiment 
was taken with a large organ (e.g. skeleton system, include 
marrow, pelvis, cartilage etc.). With the same enlargement 
size of 4.0 x 4.0 x 4.0 and a proportional weight of 0.5 for 
the proportional and the hybrid method, modeling parame- 
ters between the three samplings are given in Table 2 and the 
visualization comparison is illustrated in Fig. 4. 

In the hybrid sampling, the enlargement size of 4.0 x 4.0 x 
4.0 was divided to 2.0 x 2.0 x 2.0 for both the nearest and 
weighted. The time expenditure by the hybrid sampling was 
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Fig. 3. (Color online) Heart and face model sampling comparisons. (a) Heart model by the nearest sampling. (b) Heart model by the 
proportional sampling. (c) Face model by the nearest sampling. (d) Face model by the proportional sampling 


TABLE 2. The comparison of the skeleton system reconstruction 


Comparison Nearest 
Enlargement factor 4.0 4.0 4.0 
Proportional weight 

Voxel size (mm?) 0.6 x 0.6 x 2.0 
Times (s) 142 

Number of triangular faces 2504 154 
Holes 31635 

Model size (MB) 137.62 


about 18.15% of that by proportional. In addition, the hole 
number of the model produced by the hybrid was about 16.43 
percent of that produced by the nearest. It is advantageous to 
create the model accurately and quickly for visualization [11, 
12] .The visualization of the model constructed by the three 
sampling is shown in Fig. 4. 


Fig. 4. (Color online) Visualization comparisons between three sam- 
pling methods. (a) Skeleton system of upper half of the body by the 
nearest sampling. (b) Skeleton system of upper half of the body by 
the proportional sampling. (c) Skeleton system of upper half of the 
body by the hybrid sampling. 


The nearest sampling may cause severe distortion 
(Fig. 4(a)) when the sampling enlargement factor becomes 
larger. The hybrid sampling, therefore, is an appropriate 
method to make up for such deficiency and enhance the pro- 
cess. It is a flexible sampling method to be evolved from the 
one to another with a distributed enlargement factor. The 
combined criterion depends on the sampling ratio, dataset 
size, computer performance etc. Generally, considering the 
visual impact, the enlargement size allocated for the nearest 
sampling shall not be larger than 2.0 x 2.0 x 2.0. This can 
largely reduce the data size for fast processing with current 
configuration of a personal computer. For the demanding of 


Weighted Hybrid 

4.0 4.0 4.0 4.0 4.0 4.0 

0.5 0.5 

0.6 x 0.6 x 2.0 0.6 x 0.6 x 2.0 
2374 431 

2512610 2623 755 

4245 5200 

124.74 130.39 


more accuracy of the model (Figs. 4(b) and 4(c)), the enlarge- 
ment factor allocated for the weighted sampling should be 
larger. Because there are many organs contained in the whole 
phantom, the accepted time to visualize and check the geom- 
etry of the phantom can be satisfied by the hybrid sampling. 

From Table 2, the percent of holes to faces of the near- 
est, the weighted and the hybrid is 1.26%, 0.17% and 0.20%, 
respectively. Compared with the visualization of Fig. 4, the 
percent of about 0.20% with the model can be as a reference 
acceptance criteria for better visualization. 

By the hybrid sampling, the voxels are enlarged by a factor 
of 4.0 x 4.0 x 4.0, the whole skeleton system can be shown 
in Fig. 5. 


C. Application in Monte Carlo computational geometries 


For dose calculation in lattice geometry, the methods pre- 
sented in this paper can be applied [13] . In this case for 
Monte Carlo simulation, the pixel matrix of the whole phan- 
tom should be generated. 

The weights definition can be designed by the material, 
density or other biological property. Next, accumulate the 
weights from the resultant coordinate matrix (e.g. y), locate 
the voxel in the pixel matrix satisfying the weighted factor, 
and assign it with the identification number. Finally, the pixel 
matrixes of body can be transferred to lattice representation 
geometry in MCNP codes [14, 15]. 

However, in addition to the lattice geometry [16, 17], the 
polygon-surface model, which is constructed by Marching 
Cubes [18], can hardly avoid surface intersections and holes 
in quality. For dose calculations with deformable phan- 
toms [19, 20], repairing/improvement (e.g. filling, relax- 
ation) is needed, so as to confirm no imperfection or interfer- 
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Fig. 5. (Color online) High-accuracy skeleton system by the hybrid 
sampling. 
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IV. CONCLUSION 


In this paper, a hybrid voxel sampling algorithm, which 
merged the nearest sampling and a weighted sampling, was 
put forward. Through discussions of the model construction, 
the three sampling methods are able to deal with mass voxels. 
For accurate sampling, the weighted and the hybrid sampling 
are effective. As a flexible sampling, the hybrid sampling is 
able to show high-performance especially on construction of 
large organs. It is a valid sampling method on the visualiza- 
tion of the (MC) geometries of human voxel phantom. 

Furthermore, based on these methods, the geometries for 
Monte Carlo dose calculation were then analyzed. The phan- 
tom of polygon-surface geometry or lattice representation ge- 
ometry is able to be transferred to MC codes for dose cal- 
culations. Finally, by the proposed sampling, more accurate 
CAD (Computer Aided Design) (STEP etc.) model can be 
constructed in the field of medical diagnosis and other scien- 
tific researches. 
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